% I think the results in the paper are based on
% "script_ratio_rand_exp_fun", which does NOT average across "other" good
% j, whereas results here use "script_ratio_rand_exp_fun_average", where we
% do average over j. They are very close anyways.

clear all
clc

% num_cores = str2num(getenv('SLURM_CPUS_PER_TASK'));
% c = parcluster;
% poolobj = parpool(c, num_cores);

ratio_trimmed_mean_all = [];
trimmed_mean_ratio_all = [];
mean_ratio_all = [];
differ_all = [];

theta_NPD_bs_all = [];

for j=1:5
    
   tt = num2str(j); 
   temp = strcat('m3_4_0417_bs_',tt);
   temp2 = strcat(temp,'.mat');
   load(temp2)
   
   ratio_trimmed_mean_all = [ratio_trimmed_mean_all;
                             ratio_trimmed_mean_bs];

   trimmed_mean_ratio_all = [trimmed_mean_ratio_all;
                             trimmed_mean_ratio_bs];

   differ_all = [differ_all; ratio_trimmed_mean_bs-trimmed_mean_ratio_bs];

   mean_ratio_all = [mean_ratio_all;
                     mean_ratio_bs];
    
   theta_NPD_bs_all = [theta_NPD_bs_all;
                       theta_NPD_bs];
                         
end

nanstd(ratio_trimmed_mean_all)
nanstd(trimmed_mean_ratio_all)
nanstd(mean_ratio_all)
nanstd(differ_all)

%% Get bootstrap indices

indices = NaN(N,1000);
N = size(data_wide,1);
Nbs = 200;

counter = 1;

for i=1:5
    for bs=1:Nbs
    
    if i==1
        aa=1;
    else
        aa = i*1000;
    end

    rng(aa*bs)
    indices(:,counter) = randsample(N,N,true);

    counter = counter+1;

    end
end

writematrix(indices,'indices_bootstrap.csv')


load('data.mat')

% data = data(data.IsRevealed==0,:);

J = 3;    % # of goods

id_unique = unique(data.choice);
ngoods = zeros(length(id_unique),1);

for i=1:length(id_unique)
    ngoods(i) = sum(data.choice==id_unique(i));
end

id_unique = id_unique(ngoods==3);

data = data(ismember(data.choice,id_unique),:);

N = length(id_unique);

data_wide.choice = zeros(N,1);
data_wide.x = zeros(N,J);
data_wide.z = zeros(N,J);

for i=1:N
    ind = (data.choice==id_unique(i));
    data_wide.choice(i) = find(data.chosen(ind));
    data_wide.x(i,:) = data.price(ind)';
    data_wide.z(i,:) = data.discount(ind)';
end

data_wide = struct2table(data_wide);

data_wide.x = -data_wide.x;
data_wide.z = -data_wide.z;

[data_trans] = transform_data_01_exp(data_wide,J);

m_regr = [3 4 3*ones(1,J-2) 4 4 3*ones(1,J-2)];

[combos_all_unique] = combinations_10(m_regr,J);

% % constraints
% 
% [Aineq,bineq] = constraint_mat_theta_2(combos_all_unique,J);
% 
% Aeq = [];
% beq = [];



%% get estimate

load('m3_4.mat')

ntheta = size(theta_NPD,1);

ratio_trimmed_mean = script_ratio_rand_exp_fun_average(theta_NPD,m_regr,data_wide,J,ntheta,combos_all_unique)

%% get s.e.

Nbs = size(theta_NPD_bs_all,1);

for bs=1:Nbs
        
    ratio_trimmed_mean_bs(bs) = script_ratio_rand_exp_fun_average(theta_NPD_bs_all(bs,:)',m_regr,data_wide,J,ntheta,combos_all_unique);
    bs
end

save('m3_results.mat','ratio_trimmed_mean','ratio_trimmed_mean_bs')


